Celem analizy było określenie głównych przyczyn stopniowego zmniejszania się długości śledzi oceanicznych wyławianych w Europie. Analiza doprowadziła do wniosków, że istotny wpływ na długość miała temperatura przy powierzchni wody i to ona mogła doprowadzić do zmian długości. Duży wpływ na długość śledzia miały również wartości związane z dostępnością planktonu, co jest logiczne gdyż plankton jest pożywieniem dla śledzi. Wszystkie te wartości mogą jednak być związane z oscylacją północnoatlantycką, która jest zjawiskiem meteorologicznym i może mieć związek z temperaturą wody przy powierzchni i poprzez cyrkulację wody oceanicznej z dostępnością planktonu. Wartość ta jest przeciwnie skorelowana z długością, jednakże żaden regresor nie uznał jej za istotną. W danych brakuje informacji o roku pomiaru co utrudnia analizę, jednakże można ją wyłuskać z atrybutów rocznego narybku i rocznego natężenia połowów w rejonie.
Celem projektu jest określenie głównych przyczyn stopniowego zmniejszania się długości śledzi oceanicznych wyławionych w Europie.
Zbiór danych do analizy składa się z pomiarów śledzi i warunków w jakich żyją z ostatnich 60 lat. Dane były pobierane z połowów komercyjnych jednostek. W ramach połowu jednej jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi.
Wszystkie dane, są danymi liczbowymi, nie występują żadne dane kategoryczne.
library(dplyr)
library(ggplot2)
library(GGally) # ggpairplot
library(DT)
library(mice) # for missing data
library(ggpubr) # ggarange
library(corrplot) # ggarange
library(plotly)
library(caret)
library(gbm)
library(doMC)
columns_names <- c("length" = "długość złowionego śledzia [cm]",
"cfin1" = "dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1]",
"cfin2" = "dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2]",
"chel1" = "dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1]",
"chel2" = "dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2]",
"lcop1" = "dostępność planktonu [zagęszczenie widłonogów gat. 1]",
"lcop2" = "dostępność planktonu [zagęszczenie widłonogów gat. 2]",
"fbar" = "natężenie połowów w regionie [ułamek pozostawionego narybku]",
"recr" = "roczny narybek [liczba śledzi]",
"cumf" = "łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku]",
"totaln" = "łączna liczba ryb złowionych w ramach połowu [liczba śledzi]",
"sst" = "temperatura przy powierzchni wody [°C]",
"sal" = "poziom zasolenia wody [Knudsen ppt]",
"xmonth" = "miesiąc połowu [numer miesiąca]",
"nao" = "oscylacja północnoatlantycka [mb]" )
prettyTable <- function(table_df, round_columns=numeric(), round_digits=2) {
DT::datatable(table_df, style="bootstrap", filter = "top", rownames = FALSE, extensions = "Buttons", options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print'))) %>%
formatRound(round_columns, round_digits)
}
Pierwszym krokiem jest wczytanie danych oraz bliższe przyjżenie się im. Ponieważ w pliku csv brakujące rekordy zapisane są jako “?” trzeba je poprawnie wczytać.Brakujące dane występują w kolumnach cfin1, cfin2, chel1, chel2, lcop1, lcop2 i lsst. Jako że kolumna x to po prostu numer pomiaru, nie będzie nam potrzebna.
set.seed(23)
data <- read.csv("sledzie.csv", na.strings = c("?"))
lapply(data, function(x) any(is.na(x)))
## $X
## [1] FALSE
##
## $length
## [1] FALSE
##
## $cfin1
## [1] TRUE
##
## $cfin2
## [1] TRUE
##
## $chel1
## [1] TRUE
##
## $chel2
## [1] TRUE
##
## $lcop1
## [1] TRUE
##
## $lcop2
## [1] TRUE
##
## $fbar
## [1] FALSE
##
## $recr
## [1] FALSE
##
## $cumf
## [1] FALSE
##
## $totaln
## [1] FALSE
##
## $sst
## [1] TRUE
##
## $sal
## [1] FALSE
##
## $xmonth
## [1] FALSE
##
## $nao
## [1] FALSE
data <- subset(data, select = -c(X))
summary(data)
## length cfin1 cfin2 chel1
## Min. :19.0 Min. : 0.0000 Min. : 0.0000 Min. : 0.000
## 1st Qu.:24.0 1st Qu.: 0.0000 1st Qu.: 0.2778 1st Qu.: 2.469
## Median :25.5 Median : 0.1111 Median : 0.7012 Median : 5.750
## Mean :25.3 Mean : 0.4458 Mean : 2.0248 Mean :10.006
## 3rd Qu.:26.5 3rd Qu.: 0.3333 3rd Qu.: 1.7936 3rd Qu.:11.500
## Max. :32.5 Max. :37.6667 Max. :19.3958 Max. :75.000
## NA's :1581 NA's :1536 NA's :1555
## chel2 lcop1 lcop2 fbar
## Min. : 5.238 Min. : 0.3074 Min. : 7.849 Min. :0.0680
## 1st Qu.:13.427 1st Qu.: 2.5479 1st Qu.:17.808 1st Qu.:0.2270
## Median :21.673 Median : 7.0000 Median :24.859 Median :0.3320
## Mean :21.221 Mean : 12.8108 Mean :28.419 Mean :0.3304
## 3rd Qu.:27.193 3rd Qu.: 21.2315 3rd Qu.:37.232 3rd Qu.:0.4560
## Max. :57.706 Max. :115.5833 Max. :68.736 Max. :0.8490
## NA's :1556 NA's :1653 NA's :1591
## recr cumf totaln sst
## Min. : 140515 Min. :0.06833 Min. : 144137 Min. :12.77
## 1st Qu.: 360061 1st Qu.:0.14809 1st Qu.: 306068 1st Qu.:13.60
## Median : 421391 Median :0.23191 Median : 539558 Median :13.86
## Mean : 520367 Mean :0.22981 Mean : 514973 Mean :13.87
## 3rd Qu.: 724151 3rd Qu.:0.29803 3rd Qu.: 730351 3rd Qu.:14.16
## Max. :1565890 Max. :0.39801 Max. :1015595 Max. :14.73
## NA's :1584
## sal xmonth nao
## Min. :35.40 Min. : 1.000 Min. :-4.89000
## 1st Qu.:35.51 1st Qu.: 5.000 1st Qu.:-1.89000
## Median :35.51 Median : 8.000 Median : 0.20000
## Mean :35.51 Mean : 7.258 Mean :-0.09236
## 3rd Qu.:35.52 3rd Qu.: 9.000 3rd Qu.: 1.63000
## Max. :35.61 Max. :12.000 Max. : 5.08000
##
paste0("Ilość badanych przypadków: ", nrow(data), ", ilość zmiennych: ", ncol(data))
## [1] "Ilosc badanych przypadków: 52582, ilosc zmiennych: 15"
head(data, 10)
## length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr
## 1 23.0 0.02778 0.27785 2.46875 NA 2.54787 26.35881 0.356 482831
## 2 22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 3 25.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 4 25.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 5 24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 6 22.0 0.02778 0.27785 2.46875 21.43548 2.54787 NA 0.356 482831
## 7 24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 8 23.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 9 22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 10 22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## cumf totaln sst sal xmonth nao
## 1 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 2 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 3 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 4 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 5 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 6 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 7 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 8 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 9 0.3059879 267380.8 14.30693 35.51234 7 2.8
## 10 0.3059879 267380.8 14.30693 35.51234 7 2.8
W celu poradzenia sobie z brakującymi danymi zostaje wykorzystana metoda MICE (Multivariate Imputation via Chained Equations). Metoda ta pozwala na pozbycie się brakujących danych z utrzymaniem podobnego rozkład
plot_histograms <- function(data, columns=NULL, ncol=1, nrow=1, createName=function(x) x, bins=30){
if(is.null(columns)){
columns <- names(data)
}
i <- 1
plots <- list()
for (column_name in columns){
p <- ggplot(data, aes_string(x = column_name)) +
geom_histogram(color="black", fill="orange", bins=bins) +
# geom_density() +
ggtitle(createName(column_name)) +
ylab("Liczba obserwacji")
plots[[i]] <- p
i <- i + 1
}
figure <- ggarrange(plotlist = plots, ncol = ncol, nrow = nrow)
figure
}
missing_data_columns <- c("cfin1", "cfin2", "chel1", "chel2", "lcop1", "lcop2", "sst")
# plot_histograms(data, columns = missing_data_columns, createName = function(x) paste0("Histogram ", x, " przed pozbyciem się brakujących danych"), nrow=2)
# mice_plot <- aggr(data, col=c('navyblue','yellow'),
# numbers=TRUE, sortVars=TRUE,
# labels=names(data), cex.axis=.7,
# gap=3, ylab=c("Missing data","Pattern"))
completed_data <- data %>%
mice(m=3, method = 'cart', seed = 63) %>%
complete(2)
##
## iter imp variable
## 1 1 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 1 2 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 1 3 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 2 1 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 2 2 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 2 3 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 3 1 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 3 2 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 3 3 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 4 1 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 4 2 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 4 3 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 5 1 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 5 2 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
## 5 3 cfin1 cfin2 chel1 chel2 lcop1 lcop2 sst
# plot_histograms(completed_data, columns = missing_data_columns, createName = function(x) paste0("Histogram ", x, " po pozbyciu się brakujących danych"), nrow=2)
Podsumowanie danych
summary(completed_data)
## length cfin1 cfin2 chel1
## Min. :19.0 Min. : 0.0000 Min. : 0.0000 Min. : 0.000
## 1st Qu.:24.0 1st Qu.: 0.0000 1st Qu.: 0.2778 1st Qu.: 2.469
## Median :25.5 Median : 0.1111 Median : 0.7012 Median : 5.750
## Mean :25.3 Mean : 0.4460 Mean : 2.0255 Mean :10.000
## 3rd Qu.:26.5 3rd Qu.: 0.3333 3rd Qu.: 1.7936 3rd Qu.:11.500
## Max. :32.5 Max. :37.6667 Max. :19.3958 Max. :75.000
## chel2 lcop1 lcop2 fbar
## Min. : 5.238 Min. : 0.3074 Min. : 7.849 Min. :0.0680
## 1st Qu.:13.427 1st Qu.: 2.5479 1st Qu.:17.808 1st Qu.:0.2270
## Median :21.673 Median : 7.0000 Median :24.859 Median :0.3320
## Mean :21.219 Mean : 12.8060 Mean :28.422 Mean :0.3304
## 3rd Qu.:27.193 3rd Qu.: 21.2315 3rd Qu.:37.232 3rd Qu.:0.4560
## Max. :57.706 Max. :115.5833 Max. :68.736 Max. :0.8490
## recr cumf totaln sst
## Min. : 140515 Min. :0.06833 Min. : 144137 Min. :12.77
## 1st Qu.: 360061 1st Qu.:0.14809 1st Qu.: 306068 1st Qu.:13.60
## Median : 421391 Median :0.23191 Median : 539558 Median :13.86
## Mean : 520367 Mean :0.22981 Mean : 514973 Mean :13.87
## 3rd Qu.: 724151 3rd Qu.:0.29803 3rd Qu.: 730351 3rd Qu.:14.16
## Max. :1565890 Max. :0.39801 Max. :1015595 Max. :14.73
## sal xmonth nao
## Min. :35.40 Min. : 1.000 Min. :-4.89000
## 1st Qu.:35.51 1st Qu.: 5.000 1st Qu.:-1.89000
## Median :35.51 Median : 8.000 Median : 0.20000
## Mean :35.51 Mean : 7.258 Mean :-0.09236
## 3rd Qu.:35.52 3rd Qu.: 9.000 3rd Qu.: 1.63000
## Max. :35.61 Max. :12.000 Max. : 5.08000
Histogramy wartości atrybutów w zbiorze danych Jak można zauwyć wartości długości śledzia przyjmują rozkład normalny z średnią na poziomie 25.3 cm i przyjmuje wartości w granicach 19-32.5cm.
plot_histograms(completed_data, createName = function(x) paste0("Histogram - ", columns_names[x]), nrow=1, bins=40)
## $`1`
##
## $`2`
##
## $`3`
##
## $`4`
##
## $`5`
##
## $`6`
##
## $`7`
##
## $`8`
##
## $`9`
##
## $`10`
##
## $`11`
##
## $`12`
##
## $`13`
##
## $`14`
##
## $`15`
##
## attr(,"class")
## [1] "list" "ggarrange"
Poniżej została przedstawiona korelacja pomiędzy atrybutami. Można zauważyć, że występuje wysoka korelacja pomiędzy parami atrybutów chel1-lcop1, chel2-lcop2 oraz fbar-cumf. Korelacja może występować, ponieważ, dwie pierwsze pary związane są z dostępnością zaś ostatnia para korelacji może występować, ponieważ obydwa atrybuty związane są z natężeniem popłowu w regionie. Można te pary zredukować, aby przeciwdziałać klątwie wielowymiarowości - jednakże tutaj zostaną one zostawione bez zmian.
corrplot(cor(completed_data), type = "upper", tl.col = "black", tl.srt = 45, method="number")
ggplot(completed_data, aes(x=chel1, y=lcop1)) +
geom_point() +
geom_smooth(method=lm) +
ggtitle("Zależność pomiędzy atrybutem chel1 i lcop1")
linear_model_1 <- lm(chel1 ~ lcop1, data = completed_data)
ggplot(completed_data, aes(x=chel2, y=lcop2)) +
geom_point() +
geom_smooth(method=lm) +
ggtitle("Zależność pomiędzy atrybutem chel2 i lcop2")
linear_model_2 <- lm(chel2 ~ lcop2, data = completed_data)
ggplot(completed_data, aes(x=fbar, y=cumf)) +
geom_point() +
geom_smooth(method=lm) +
ggtitle("Zależność pomiędzy atrybutem fbar i cumf")
linear_model_3 <- lm(fbar ~ cumf, data = completed_data)
Poniżej przedstawiona została analiza długości śledzia na przestrzeni czasu. Jednakże czas nie był przedstawiony jawnie, ale wiadomo, że dane są ustawione chronologicznie dzięki czemu, zbiór można podzielić i na każdym z podziałów policzyć średnią. Można zauważyć, że na początku długość śledzia na przestrzeni czasu wzrastała od średniej wielkości w oklicach 24.5 cm do wielkości wielkości przekraczającej 27cm. Jednakże później trend się odwrócił i ostatecznie na sam koniec średnia ta była poniżej 24 cm.
aggregate_data <- function(data, cases_per_aggregate){
new_aggregated_data <- aggregate(data,
list(rep(1:(nrow(data)%/%cases_per_aggregate+1), each=cases_per_aggregate, len=nrow(data))),
mean) # Try using mean without outliners
new_aggregated_data$time_period = c(1:nrow(new_aggregated_data))
new_aggregated_data
}
aggregated_data <- aggregate_data(completed_data, 200)
ggplot(aggregated_data, aes(time_period, length)) +
geom_point() +
geom_smooth() +
ggtitle("Wykres zmiany długości śledzia na przestrzeni czasu") +
ylab("Długość śledzia") +
xlab("Przedział czasowy")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
plot_ly(data = aggregated_data, x = ~time_period, y = ~length,
text = ~paste("Length: ", length))
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
# %>% layout(
# scene = list(
# xaxis = list(title = "Przedział czasowy"),
# yaxis = list(title = "Długość śledzia")
# )
# )
train_index <- createDataPartition(completed_data$length, p = .7,
list = FALSE,
times = 1)
training <- completed_data[ train_index,]
testing <- completed_data[-train_index,]
r_squared <- function (x, y) cor(x, y) ^ 2
registerDoMC(cores=8)
ctrl <- trainControl(
method = "repeatedcv",
number = 5,
repeats = 3
)
evaluate_model <- function(model, test_data) {
print(model)
plot <- ggplot(varImp(model))
print(varImp(model))
print(plot)
predicted <- predict(model, newdata = test_data)
rmse_ <- RMSE(predicted, test_data$length)
rsquared <- r_squared(predicted, test_data$length)
print(paste0("RMSE na zbiorze testowym: ", rmse_, " R-Squared na zbiorze testowym ", rsquared))
c(rmse_, rsquared)
}
modelLookup(model='glmnet')
## model parameter label forReg forClass probModel
## 1 glmnet alpha Mixing Percentage TRUE TRUE TRUE
## 2 glmnet lambda Regularization Parameter TRUE TRUE TRUE
glmnet_grid = expand.grid(
alpha = 0:1,
lambda = seq(0.0001, 1, length = 100)
)
set.seed(23)
glmnet <- train(length ~ .,
data = training,
method = "glmnet",
preProc = c("center", "scale"),
metric = "RMSE",
trControl = ctrl,
tuneGrid = glmnet_grid,
importance=TRUE)
## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info =
## trainInfo, : There were missing values in resampled performance measures.
glmnet_metrics = evaluate_model(glmnet, testing)
## glmnet
##
## 36810 samples
## 14 predictor
##
## Pre-processing: centered (14), scaled (14)
## Resampling: Cross-Validated (5 fold, repeated 3 times)
## Summary of sample sizes: 29447, 29449, 29449, 29447, 29448, 29448, ...
## Resampling results across tuning parameters:
##
## alpha lambda RMSE Rsquared MAE
## 0 0.0001 1.370182 0.3143406 1.091332
## 0 0.0102 1.370182 0.3143406 1.091332
## 0 0.0203 1.370182 0.3143406 1.091332
## 0 0.0304 1.370182 0.3143406 1.091332
## 0 0.0405 1.370182 0.3143406 1.091332
## 0 0.0506 1.370182 0.3143406 1.091332
## 0 0.0607 1.370182 0.3143406 1.091332
## 0 0.0708 1.370182 0.3143406 1.091332
## 0 0.0809 1.370975 0.3137485 1.091997
## 0 0.0910 1.372204 0.3128238 1.093068
## 0 0.1011 1.373432 0.3118982 1.094216
## 0 0.1112 1.374655 0.3109730 1.095383
## 0 0.1213 1.375870 0.3100524 1.096521
## 0 0.1314 1.377076 0.3091374 1.097612
## 0 0.1415 1.378265 0.3082344 1.098650
## 0 0.1516 1.379429 0.3073522 1.099638
## 0 0.1617 1.380579 0.3064784 1.100609
## 0 0.1718 1.381721 0.3056095 1.101582
## 0 0.1819 1.382818 0.3047808 1.102528
## 0 0.1920 1.383915 0.3039490 1.103488
## 0 0.2021 1.384983 0.3031414 1.104450
## 0 0.2122 1.386037 0.3023462 1.105412
## 0 0.2223 1.387075 0.3015634 1.106367
## 0 0.2324 1.388090 0.3008004 1.107301
## 0 0.2425 1.389083 0.3000564 1.108214
## 0 0.2526 1.390069 0.2993171 1.109116
## 0 0.2627 1.391016 0.2986134 1.109975
## 0 0.2728 1.391980 0.2978926 1.110840
## 0 0.2829 1.392891 0.2972216 1.111651
## 0 0.2930 1.393809 0.2965424 1.112462
## 0 0.3031 1.394718 0.2958696 1.113262
## 0 0.3132 1.395580 0.2952419 1.114014
## 0 0.3233 1.396456 0.2945998 1.114773
## 0 0.3334 1.397314 0.2939733 1.115509
## 0 0.3435 1.398135 0.2933822 1.116206
## 0 0.3536 1.398968 0.2927788 1.116905
## 0 0.3637 1.399796 0.2921783 1.117598
## 0 0.3738 1.400574 0.2916267 1.118248
## 0 0.3839 1.401362 0.2910647 1.118907
## 0 0.3940 1.402160 0.2904927 1.119576
## 0 0.4041 1.402916 0.2899607 1.120212
## 0 0.4142 1.403659 0.2894412 1.120834
## 0 0.4243 1.404410 0.2889128 1.121458
## 0 0.4344 1.405166 0.2883789 1.122081
## 0 0.4445 1.405873 0.2878921 1.122659
## 0 0.4546 1.406576 0.2874098 1.123228
## 0 0.4647 1.407286 0.2869200 1.123797
## 0 0.4748 1.408002 0.2864233 1.124366
## 0 0.4849 1.408681 0.2859627 1.124903
## 0 0.4950 1.409343 0.2855180 1.125424
## 0 0.5051 1.410012 0.2850668 1.125947
## 0 0.5152 1.410687 0.2846092 1.126474
## 0 0.5253 1.411360 0.2841532 1.127002
## 0 0.5354 1.411988 0.2837410 1.127506
## 0 0.5455 1.412617 0.2833277 1.128024
## 0 0.5556 1.413252 0.2829089 1.128557
## 0 0.5657 1.413892 0.2824845 1.129103
## 0 0.5758 1.414529 0.2820619 1.129652
## 0 0.5859 1.415124 0.2816807 1.130171
## 0 0.5960 1.415716 0.2813025 1.130691
## 0 0.6061 1.416312 0.2809196 1.131218
## 0 0.6162 1.416913 0.2805321 1.131752
## 0 0.6263 1.417518 0.2801402 1.132292
## 0 0.6364 1.418103 0.2797678 1.132816
## 0 0.6465 1.418660 0.2794233 1.133318
## 0 0.6566 1.419220 0.2790749 1.133823
## 0 0.6667 1.419784 0.2787226 1.134330
## 0 0.6768 1.420352 0.2783663 1.134839
## 0 0.6869 1.420923 0.2780063 1.135349
## 0 0.6970 1.421481 0.2776595 1.135845
## 0 0.7071 1.422006 0.2773437 1.136314
## 0 0.7172 1.422535 0.2770248 1.136783
## 0 0.7273 1.423068 0.2767024 1.137252
## 0 0.7374 1.423603 0.2763767 1.137722
## 0 0.7475 1.424142 0.2760476 1.138192
## 0 0.7576 1.424682 0.2757175 1.138660
## 0 0.7677 1.425196 0.2754119 1.139106
## 0 0.7778 1.425691 0.2751242 1.139535
## 0 0.7879 1.426190 0.2748338 1.139964
## 0 0.7980 1.426691 0.2745405 1.140394
## 0 0.8081 1.427195 0.2742443 1.140825
## 0 0.8182 1.427702 0.2739453 1.141257
## 0 0.8283 1.428211 0.2736438 1.141689
## 0 0.8384 1.428708 0.2733533 1.142112
## 0 0.8485 1.429176 0.2730912 1.142511
## 0 0.8586 1.429646 0.2728278 1.142911
## 0 0.8687 1.430117 0.2725620 1.143313
## 0 0.8788 1.430591 0.2722939 1.143717
## 0 0.8889 1.431068 0.2720233 1.144122
## 0 0.8990 1.431546 0.2717502 1.144528
## 0 0.9091 1.432027 0.2714752 1.144935
## 0 0.9192 1.432498 0.2712089 1.145333
## 0 0.9293 1.432941 0.2709697 1.145707
## 0 0.9394 1.433383 0.2707307 1.146079
## 0 0.9495 1.433828 0.2704896 1.146452
## 0 0.9596 1.434274 0.2702466 1.146825
## 0 0.9697 1.434723 0.2700015 1.147198
## 0 0.9798 1.435173 0.2697543 1.147572
## 0 0.9899 1.435626 0.2695050 1.147946
## 0 1.0000 1.436079 0.2692550 1.148319
## 1 0.0001 1.363552 0.3192588 1.084409
## 1 0.0102 1.365472 0.3177033 1.085918
## 1 0.0203 1.369221 0.3148282 1.089374
## 1 0.0304 1.374165 0.3110152 1.093505
## 1 0.0405 1.379621 0.3068474 1.099438
## 1 0.0506 1.386537 0.3011004 1.106289
## 1 0.0607 1.392899 0.2959747 1.112238
## 1 0.0708 1.399495 0.2905558 1.118008
## 1 0.0809 1.407027 0.2838729 1.124723
## 1 0.0910 1.415296 0.2760606 1.132537
## 1 0.1011 1.423924 0.2675317 1.140801
## 1 0.1112 1.433152 0.2578987 1.149419
## 1 0.1213 1.441206 0.2494461 1.156718
## 1 0.1314 1.446138 0.2449423 1.161262
## 1 0.1415 1.448712 0.2436123 1.163650
## 1 0.1516 1.451009 0.2427450 1.165754
## 1 0.1617 1.453433 0.2417993 1.167946
## 1 0.1718 1.456002 0.2407395 1.170163
## 1 0.1819 1.458721 0.2395443 1.172424
## 1 0.1920 1.461590 0.2382005 1.174756
## 1 0.2021 1.464608 0.2366931 1.177180
## 1 0.2122 1.467774 0.2350070 1.180064
## 1 0.2223 1.471086 0.2331238 1.183065
## 1 0.2324 1.474546 0.2310234 1.186110
## 1 0.2425 1.478145 0.2286967 1.189240
## 1 0.2526 1.481634 0.2265382 1.192206
## 1 0.2627 1.485212 0.2242264 1.195414
## 1 0.2728 1.488566 0.2223093 1.198411
## 1 0.2829 1.491748 0.2207225 1.201208
## 1 0.2930 1.495038 0.2189623 1.204101
## 1 0.3031 1.498436 0.2170113 1.207070
## 1 0.3132 1.501941 0.2148499 1.210146
## 1 0.3233 1.505550 0.2124620 1.213486
## 1 0.3334 1.509226 0.2098939 1.216839
## 1 0.3435 1.513005 0.2070578 1.220220
## 1 0.3536 1.516780 0.2041574 1.223578
## 1 0.3637 1.519431 0.2036212 1.225840
## 1 0.3738 1.522063 0.2032341 1.228081
## 1 0.3839 1.524757 0.2028211 1.230388
## 1 0.3940 1.527385 0.2026949 1.232609
## 1 0.4041 1.530021 0.2026949 1.234872
## 1 0.4142 1.532719 0.2026949 1.237169
## 1 0.4243 1.535479 0.2026949 1.239466
## 1 0.4344 1.538300 0.2026949 1.241771
## 1 0.4445 1.541182 0.2026949 1.244114
## 1 0.4546 1.544125 0.2026949 1.246473
## 1 0.4647 1.547128 0.2026949 1.248842
## 1 0.4748 1.550191 0.2026949 1.251314
## 1 0.4849 1.553314 0.2026949 1.253911
## 1 0.4950 1.556496 0.2026949 1.256537
## 1 0.5051 1.559737 0.2026949 1.259215
## 1 0.5152 1.563036 0.2026949 1.262035
## 1 0.5253 1.566394 0.2026949 1.264921
## 1 0.5354 1.569809 0.2026949 1.267861
## 1 0.5455 1.573282 0.2026949 1.270836
## 1 0.5556 1.576812 0.2026949 1.273878
## 1 0.5657 1.580399 0.2026949 1.276973
## 1 0.5758 1.584042 0.2026949 1.280119
## 1 0.5859 1.587741 0.2026949 1.283340
## 1 0.5960 1.591495 0.2026949 1.286631
## 1 0.6061 1.595305 0.2026949 1.290000
## 1 0.6162 1.599169 0.2026949 1.293437
## 1 0.6263 1.603087 0.2026949 1.296960
## 1 0.6364 1.607060 0.2026949 1.300528
## 1 0.6465 1.611086 0.2026949 1.304142
## 1 0.6566 1.615165 0.2026949 1.307781
## 1 0.6667 1.619297 0.2026949 1.311441
## 1 0.6768 1.623481 0.2026949 1.315195
## 1 0.6869 1.627717 0.2026949 1.318972
## 1 0.6970 1.632005 0.2026949 1.322749
## 1 0.7071 1.636343 0.2026949 1.326527
## 1 0.7172 1.640733 0.2026949 1.330304
## 1 0.7273 1.645172 0.2026949 1.334081
## 1 0.7374 1.649614 0.2015148 1.337818
## 1 0.7475 1.652483 0.1944352 1.340208
## 1 0.7576 1.652485 NaN 1.340209
## 1 0.7677 1.652485 NaN 1.340209
## 1 0.7778 1.652485 NaN 1.340209
## 1 0.7879 1.652485 NaN 1.340209
## 1 0.7980 1.652485 NaN 1.340209
## 1 0.8081 1.652485 NaN 1.340209
## 1 0.8182 1.652485 NaN 1.340209
## 1 0.8283 1.652485 NaN 1.340209
## 1 0.8384 1.652485 NaN 1.340209
## 1 0.8485 1.652485 NaN 1.340209
## 1 0.8586 1.652485 NaN 1.340209
## 1 0.8687 1.652485 NaN 1.340209
## 1 0.8788 1.652485 NaN 1.340209
## 1 0.8889 1.652485 NaN 1.340209
## 1 0.8990 1.652485 NaN 1.340209
## 1 0.9091 1.652485 NaN 1.340209
## 1 0.9192 1.652485 NaN 1.340209
## 1 0.9293 1.652485 NaN 1.340209
## 1 0.9394 1.652485 NaN 1.340209
## 1 0.9495 1.652485 NaN 1.340209
## 1 0.9596 1.652485 NaN 1.340209
## 1 0.9697 1.652485 NaN 1.340209
## 1 0.9798 1.652485 NaN 1.340209
## 1 0.9899 1.652485 NaN 1.340209
## 1 1.0000 1.652485 NaN 1.340209
##
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 1 and lambda = 1e-04.
## glmnet variable importance
##
## Overall
## fbar 100.000
## cumf 95.433
## sst 54.392
## lcop1 30.601
## chel1 12.388
## totaln 11.977
## lcop2 11.788
## cfin1 11.436
## recr 9.039
## cfin2 6.426
## nao 4.945
## chel2 1.719
## xmonth 1.099
## sal 0.000
## [1] "RMSE na zbiorze testowym: 1.35944937812852 R-Squared na zbiorze testowym 0.324293505084823"
set.seed(23)
modelLookup(model='rf')
## model parameter label forReg forClass probModel
## 1 rf mtry #Randomly Selected Predictors TRUE TRUE TRUE
rf_grid <- expand.grid(mtry=c(2, 3, 4, 5, 6, 7, 8, 9, 11, 14))
random_forest <- train(length ~ .,
data = training,
method = "rf",
preProc = c("center", "scale"),
metric = "RMSE",
trControl = ctrl,
tuneGrid = rf_grid,
ntree = 10,
importance=TRUE)
rf_metrics = evaluate_model(random_forest, testing)
## Random Forest
##
## 36810 samples
## 14 predictor
##
## Pre-processing: centered (14), scaled (14)
## Resampling: Cross-Validated (5 fold, repeated 3 times)
## Summary of sample sizes: 29447, 29449, 29449, 29447, 29448, 29448, ...
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 2 1.168083 0.5005193 0.9238434
## 3 1.157629 0.5093942 0.9146492
## 4 1.151832 0.5142911 0.9088003
## 5 1.147856 0.5176484 0.9049700
## 6 1.146889 0.5185441 0.9034054
## 7 1.146685 0.5188028 0.9028100
## 8 1.148130 0.5177382 0.9036876
## 9 1.148550 0.5174491 0.9039509
## 11 1.150159 0.5162438 0.9051027
## 14 1.150568 0.5159722 0.9052872
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 7.
## rf variable importance
##
## Overall
## xmonth 100.0000
## chel1 32.0915
## cfin2 23.6747
## chel2 17.8241
## sst 12.8506
## lcop2 12.6688
## totaln 11.1845
## lcop1 9.1186
## cfin1 8.8539
## cumf 7.9361
## fbar 3.9435
## recr 2.4297
## sal 0.2044
## nao 0.0000
## [1] "RMSE na zbiorze testowym: 1.1447769654101 R-Squared na zbiorze testowym 0.520876878600763"
Podane parametry zostały znaliezone wcześniej wykorzystując grid search, jednakż trwało to bardzo długo i żeby nie ponawiać tych obliczeń parametry zostały tu wpisane na sztywno.
modelLookup(model='gbm')
## model parameter label forReg forClass
## 1 gbm n.trees # Boosting Iterations TRUE TRUE
## 2 gbm interaction.depth Max Tree Depth TRUE TRUE
## 3 gbm shrinkage Shrinkage TRUE TRUE
## 4 gbm n.minobsinnode Min. Terminal Node Size TRUE TRUE
## probModel
## 1 TRUE
## 2 TRUE
## 3 TRUE
## 4 TRUE
gbm_grid <- expand.grid(
n.trees=c(10,20,50,100,500,1000),
shrinkage=c(0.01,0.05,0.1,0.5),
n.minobsinnode = c(3,5,10),
interaction.depth=c(1,5,10)
)
set.seed(23)
params <- data.frame( shrinkage = 0.1, n.trees = 1000, interaction.depth = 10, n.minobsinnode = 10)
gbm <- train(length ~ .,
data = training,
method = "gbm",
metric = "RMSE",
tuneGrid = params,
# trControl = ctrl,
# tuneGrid = gbm_grid,
preProc = c("center", "scale"),
verbose = FALSE)
gbm_metrics = evaluate_model(gbm, testing)
## Stochastic Gradient Boosting
##
## 36810 samples
## 14 predictor
##
## Pre-processing: centered (14), scaled (14)
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 36810, 36810, 36810, 36810, 36810, 36810, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 1.148978 0.5181491 0.9045936
##
## Tuning parameter 'n.trees' was held constant at a value of 1000
## 10
## Tuning parameter 'shrinkage' was held constant at a value of
## 0.1
## Tuning parameter 'n.minobsinnode' was held constant at a value of 10
## gbm variable importance
##
## Overall
## sst 100.0000
## recr 27.6999
## xmonth 19.9096
## lcop1 14.6281
## totaln 8.4273
## cfin2 7.0833
## lcop2 6.6431
## chel2 4.8450
## cfin1 3.0231
## sal 0.8983
## nao 0.7995
## cumf 0.6333
## fbar 0.5033
## chel1 0.0000
## [1] "RMSE na zbiorze testowym: 1.13672654284413 R-Squared na zbiorze testowym 0.527527267988041"
ggplot(varImp(gbm))
Jak można zauważyć poniżej najlepszymi modelami okazały się modele oparte na algorytmach Gradient Boosting i Random Forest. Jest nieznaczna różnica pomiędzy tymi algorytmami. Nieco gorzej sprawdza się algorytm random forest. Najgorzej i poniżej oczekiwań działa algorytm Elastic Net, co może być związane z wykorzystaniem regresji liniowej wewnątrz algorytmu.
metrics <- data.frame(
regressor = c("Elastic Net","Random Forest","Gradient Boosting Regressor"),
rmse = c(glmnet_metrics[1], rf_metrics[1], gbm_metrics[1]),
rsquared = c(glmnet_metrics[2], rf_metrics[2], gbm_metrics[2])
)
metrics
## regressor rmse rsquared
## 1 Elastic Net 1.359449 0.3242935
## 2 Random Forest 1.144777 0.5208769
## 3 Gradient Boosting Regressor 1.136727 0.5275273
Analizując istotność danych dla Elastic Net trzema najistotniejszymi atrybutami są fbar, cumf i sst, a najmniej istotnymi chel2, xmonth i sal. Dla algorytmu random forest najistotniejszymi były xmonth, chel1 i cfin2, a najmniej istotne cumf, sal i mao. Najważniejszymi atrybutami dla algorytmu Gradient Boosting były sst, recr i xmonth, a najmniej istotnymi chel1, fbar i sal (wszystkie bardzo niskie).
Można się spodziewać, że atrybut xmonth jest skorelowany z długością (jest to najistotniejszy atrybut dla Random Forest i trzeci najbardziej istotny dla Gradient Boosting). Jednakże nie oznacza to, że jest to wynik wieloletniego trendu, a może na przykład wynikać z wahań związanych z okresami godowymi lub okresami zwiększonego połowu.
Można zauważyć dużą istotność atrybutu sst - najistotniejszy atrybut dla Gradient Boosting, trzeci najbardziej istotny dla Elastic Net. Na wykresie poniżej można zauważyć że średnia długość śledzia jest przeciwnie skorelowana z temperaturą wody przy powierzchni. Może to świadczyć, że jest to powód zmian długości śledzi na przestrzeni lat.
ggplot(aggregated_data, aes(x = time_period)) +
geom_point(aes(y = sst * 2, colour = "sst")) +
geom_point(aes(y = length, colour = "length")) +
scale_y_continuous(sec.axis = sec_axis(~.*0.5, name = columns_names["lcop1"])) +
scale_colour_manual(values = c("blue", "red")) +
labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
theme(legend.position = c(0.25, 0.1)) +
ggtitle("Długość śledzia wraz z temperaturą pod powierzchnią wody na przestrzeni lat")
Mimo występowania korelacji pomiędzy długością, a oscylacją północnoatlantycką, żaden z regresorów nie uznawał tego atrybutu za istotny. Mniejszą zależność widać w przypadku lcop2 (i chel2).
ggplot(aggregate_data(completed_data, 1000), aes(x = time_period)) +
geom_line(aes(y = nao + 25, colour = "nao")) +
geom_line(aes(y = length, colour = "length")) +
scale_y_continuous(sec.axis = sec_axis(~.-25, name = columns_names["nao"])) +
scale_colour_manual(values = c("blue", "red")) +
labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
theme(legend.position = c(0.25, 0.1)) +
ggtitle("Długość śledzia wraz z oscylacja północnoatlantycką na przestrzeni lat")
Często istotnym atrybutem jest atrybut lcop1 (jak i chel1 poprzez korelacje) czyli dostępność planktonu [zagęszczenie widłonogów gat. 1]. Można zauważyć podobny trend zmian tych wartości do trendu zmian wartości długości. Dodatkowo jest on skorelowany z oscylacją północnoatlantycką, przez co ten drugi atrybut może nie występować w istotnych atrybutach.
ggplot(aggregate_data(completed_data, 1000), aes(x = time_period)) +
geom_line(aes(y = lcop1 / 10 + 25 , colour = "lcop1")) +
geom_line(aes(y = length, colour = "length")) +
scale_y_continuous(sec.axis = sec_axis(~ (.- 25) * 10, name = columns_names["lcop1"])) +
scale_colour_manual(values = c("blue", "red")) +
labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
theme(legend.position = c(0.25, 0.1)) +
ggtitle("Długość śledzia wraz z zagęszczeniem widłogonów gat. 1 na przestrzeni lat")
ggplot(aggregate_data(completed_data, 1000), aes(x = time_period)) +
geom_line(aes(y = chel1 / 10 + 25 , colour = "chel1")) +
geom_line(aes(y = length, colour = "length")) +
scale_y_continuous(sec.axis = sec_axis(~ (.- 25) * 10, name = columns_names["chel1"])) +
scale_colour_manual(values = c("blue", "red")) +
labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
theme(legend.position = c(0.25, 0.1)) +
ggtitle("Długość śledzia wraz z atrubutem chel1 na przestrzeni lat")
ggplot(aggregate_data(completed_data, 1000), aes(x = time_period)) +
geom_line(aes(y = lcop2 / 10 + 22 , colour = "lcop2")) +
geom_line(aes(y = length, colour = "length")) +
scale_y_continuous(sec.axis = sec_axis(~ (.- 22) * 10, name = columns_names["lcop2"])) +
scale_colour_manual(values = c("blue", "red")) +
labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
theme(legend.position = c(0.25, 0.1)) +
ggtitle("Długość śledzia wraz z atrubutem lcop2 na przestrzeni lat")
Wartości cumf i rect wydają się nie mieć wpływu na zmianę długości śledzia na przestrzeni lat.
ggplot(aggregate_data(completed_data, 1000), aes(x = time_period)) +
geom_line(aes(y = cumf * 100 , colour = "cumf")) +
geom_line(aes(y = length, colour = "length")) +
scale_y_continuous(sec.axis = sec_axis(~. / 100, name = columns_names["cumf"])) +
scale_colour_manual(values = c("blue", "red")) +
labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
theme(legend.position = c(0.25, 0.1)) +
ggtitle("Długość śledzia wraz z rocznym natęzeniem połowów na przestrzeni lat")
ggplot(aggregate_data(completed_data, 1000), aes(x = time_period)) +
geom_line(aes(y = recr * 0.00005 , colour = "recr")) +
geom_line(aes(y = length, colour = "length")) +
scale_y_continuous(sec.axis = sec_axis(~. / 0.00005, name = columns_names["recr"])) +
scale_colour_manual(values = c("blue", "red")) +
labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
theme(legend.position = c(0.25, 0.1)) +
ggtitle("Długość śledzia wraz z rocznym narybkiem na przestrzeni lat")